#!/usr/bin/env python
# -*- coding: utf-8 -*-

###############################################################################
#                                                                             #
# RMG - Reaction Mechanism Generator                                          #
#                                                                             #
# Copyright (c) 2002-2019 Prof. William H. Green (whgreen@mit.edu),           #
# Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu)   #
#                                                                             #
# Permission is hereby granted, free of charge, to any person obtaining a     #
# copy of this software and associated documentation files (the 'Software'),  #
# to deal in the Software without restriction, including without limitation   #
# the rights to use, copy, modify, merge, publish, distribute, sublicense,    #
# and/or sell copies of the Software, and to permit persons to whom the       #
# Software is furnished to do so, subject to the following conditions:        #
#                                                                             #
# The above copyright notice and this permission notice shall be included in  #
# all copies or substantial portions of the Software.                         #
#                                                                             #
# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR  #
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,    #
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE #
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER      #
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING     #
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER         #
# DEALINGS IN THE SOFTWARE.                                                   #
#                                                                             #
###############################################################################

"""
Contains classes for working with the reaction model generated by RMG.
"""

import logging
import math
import numpy
import itertools
import gc
import os

from rmgpy.display import display
from rmgpy import settings
import rmgpy.constants as constants
from rmgpy.constraints import failsSpeciesConstraints
from rmgpy.quantity import Quantity
from rmgpy.species import Species
from rmgpy.thermo.thermoengine import submit
from rmgpy.reaction import Reaction
from rmgpy.exceptions import ForbiddenStructureException
from rmgpy.data.kinetics.depository import DepositoryReaction
from rmgpy.data.kinetics.family import KineticsFamily, TemplateReaction
from rmgpy.data.kinetics.library import KineticsLibrary, LibraryReaction

from rmgpy.kinetics import KineticsData, Arrhenius

from rmgpy.data.rmg import getDB
        
import rmgpy.data.rmg
from .react import reactAll

from pdep import PDepReaction, PDepNetwork

# generateThermoDataFromQM under the Species class imports the qm package

################################################################################

class ReactionModel:
    """
    Represent a generic reaction model. A reaction model consists of `species`,
    a list of species, and `reactions`, a list of reactions.
    """

    def __init__(self, species=None, reactions=None):
        self.species = species or []
        self.reactions = reactions or []
    
    def __reduce__(self):
        """
        A helper function used when pickling an object.
        """
        return (ReactionModel, (self.species, self.reactions))

    def merge(self, other):
        """
        Return a new :class:`ReactionModel` object that is the union of this
        model and `other`.
        """
        if not isinstance(other, ReactionModel):
            raise ValueError('Expected type ReactionModel for other parameter, got {0}'.format(other.__class__))

        # Initialize the merged model
        finalModel = ReactionModel()
        
        # Put the current model into the merged model as-is
        finalModel.species.extend(self.species)
        finalModel.reactions.extend(self.reactions)
        
        # Determine which species in other are already in self
        commonSpecies = {}; uniqueSpecies = []
        for spec in other.species:
            for spec0 in finalModel.species:
                if spec.isIsomorphic(spec0):
                    commonSpecies[spec] = spec0
                    if spec0.label not in ['Ar','N2','Ne','He']:
                        if not spec0.thermo.isIdenticalTo(spec.thermo):
                            print 'Species {0} thermo from model 1 did not match that of model 2.'.format(spec.label)
                        
                    break
            else:
                uniqueSpecies.append(spec)
        
        # Determine which reactions in other are already in self
        commonReactions = {}; uniqueReactions = []
        for rxn in other.reactions:
            for rxn0 in finalModel.reactions:
                if rxn.isIsomorphic(rxn0, eitherDirection=True):
                    commonReactions[rxn] = rxn0                    
                    if not rxn0.kinetics.isIdenticalTo(rxn.kinetics):
                        print 'Reaction {0} kinetics from model 1 did not match that of model 2.'.format(str(rxn0))
                    break
            else:
                uniqueReactions.append(rxn)
        
        # Add the unique species from other to the final model
        finalModel.species.extend(uniqueSpecies)
    
        # Renumber the unique species (to avoid name conflicts on save)
        speciesIndex = 0
        for spec in finalModel.species:
            if spec.label not in ['Ar','N2','Ne','He']:
                spec.index = speciesIndex + 1
                speciesIndex += 1
        
        # Make sure unique reactions only refer to species in the final model
        for rxn in uniqueReactions:
            for i, reactant in enumerate(rxn.reactants):
                try:
                    rxn.reactants[i] = commonSpecies[reactant]
                    if rxn.pairs:
                        for j, pair in enumerate(rxn.pairs):
                            if reactant in pair:
                                rxn.pairs[j] = (rxn.reactants[i],pair[1])
                except KeyError:
                    pass
            for i, product in enumerate(rxn.products):
                try:
                    rxn.products[i] = commonSpecies[product]
                    if rxn.pairs:
                        for j, pair in enumerate(rxn.pairs):
                            if product in pair:
                                rxn.pairs[j] = (pair[0], rxn.products[i])
                except KeyError:
                    pass
        
        # Add the unique reactions from other to the final model
        finalModel.reactions.extend(uniqueReactions)
    
        # Return the merged model
        return finalModel

################################################################################

class CoreEdgeReactionModel:
    """
    Represent a reaction model constructed using a rate-based screening
    algorithm. The species and reactions in the model itself are called the
    *core*; the species and reactions identified as candidates for inclusion in
    the model are called the *edge*. The attributes are:

    =========================  ==============================================================
    Attribute                  Description
    =========================  ==============================================================
    `core`                     The species and reactions of the current model core
    `edge`                     The species and reactions of the current model edge
    `networkDict`              A dictionary of pressure-dependent reaction networks (:class:`Network` objects) indexed by source.
    `networkList`              A list of pressure-dependent reaction networks (:class:`Network` objects)
    `networkCount`             A counter for the number of pressure-dependent networks created
    `indexSpeciesDict`         A dictionary with a unique index pointing to the species objects
    `solventName`              String describing solvent name for liquid reactions. Empty for non-liquid estimation
    =========================  ==============================================================


    """

    def __init__(self, core=None, edge=None, surface=None):
        if core is None:
            self.core = ReactionModel()
        else:
            self.core = core
        if edge is None:
            self.edge = ReactionModel()
        else:
            self.edge = edge
        if surface is None:
            self.surface = ReactionModel()
        else:
            self.surface = surface
            
        # The default tolerances mimic the original RMG behavior; no edge
        # pruning takes place, and the simulation is interrupted as soon as
        # a species flux higher than the validity
        self.networkDict = {}
        self.networkList = []
        self.networkCount = 0
        self.speciesDict = {}
        self.reactionDict = {}
        self.speciesCache = [None for i in range(4)]
        self.speciesCounter = 0
        self.reactionCounter = 0
        self.newSpeciesList = []
        self.newReactionList = []
        self.outputSpeciesList = []
        self.outputReactionList = []
        self.pressureDependence = None
        self.quantumMechanics = None
        self.verboseComments = False
        self.kineticsEstimator = 'rate rules'
        self.indexSpeciesDict = {}
        self.saveEdgeSpecies = False
        self.iterationNum = 0
        self.toleranceThermoKeepSpeciesInEdge = numpy.inf
        self.Gfmax = numpy.inf
        self.Gmax = numpy.inf
        self.Gmin = -numpy.inf
        self.minCoreSizeForPrune = 50
        self.maximumEdgeSpecies = 100000
        self.Tmax = 0
        self.reactionSystems = []
        self.newSurfaceSpcsAdd = set()
        self.newSurfaceRxnsAdd = set()
        self.newSurfaceSpcsLoss = set()
        self.newSurfaceRxnsLoss = set()
        self.solventName = ''

    def checkForExistingSpecies(self, molecule):
        """
        Check to see if an existing species contains the same
        :class:`molecule.Molecule` as `molecule`.
        Returns ``True``, `reactive`, and the matched species (if found) or
        ``False``, ``False``, and ``None`` (if not found).
        `reactive` is a boolean argument which is ``False`` if this molecule is an unrepresentative resonance structure
        of an existing species (i.e., was found to be isomorphic only by generating its unfiltered resonance structures)
        and True otherwise. It is emphasized that `reactive` relates to the :Class:`Molecule` attribute.
        """
        # Create obj to check against existing species
        # obj can be `Molecule` object or `Species` object

        # For non-cyclic molecules, obj is `Molecule` object
        # We expect it to be part of the list of isomers in a species
        # object if it has a match
        obj = molecule

        # For cyclic molecules, obj is `Species` object and aromatic resonance
        # isomers are generated.  This is due to the hysteresis of isomer generation
        # for aromatic/polyaromatic compounds: not all kekulized forms can be found
        # within the list of isomers for a species object describing a unique aromatic compound
        if molecule.isCyclic():
            obj = Species(molecule=[molecule])
            from rmgpy.molecule.resonance import generate_optimal_aromatic_resonance_structures
            aromaticIsomers = generate_optimal_aromatic_resonance_structures(molecule)
            obj.molecule.extend(aromaticIsomers)

        # First check cache and return if species is found
        for i, spec in enumerate(self.speciesCache):
            if spec is not None:
                for mol in spec.molecule:
                    if obj.isIsomorphic(mol):
                        self.speciesCache.pop(i)
                        self.speciesCache.insert(0, spec)
                        return True, True, spec

        # Return an existing species if a match is found
        formula = molecule.getFormula()
        try:
            speciesList = self.speciesDict[formula]
        except KeyError:
            return False, False, None
        for spec in speciesList:
            if spec.isIsomorphic(obj):
                self.speciesCache.pop()
                self.speciesCache.insert(0, spec)
                return True, True, spec

        # As a last resort, check using molecule.fingerprint if the object matches any existing species,
        # and if it does, generate resonance structures w/o filtration and check for isomorphism
        candidates = []
        for spec in speciesList:
            if spec.molecule[0].fingerprint == molecule.fingerprint:
                candidates.append(spec)
        if len(candidates) > 0:
            mol_copy = molecule.copy(deep=True)
            if not mol_copy.reactive:
                mol_copy.reactive = True
            structures = mol_copy.generate_resonance_structures(keep_isomorphic=False, filter_structures=False)
            for spec in candidates:
                for mol in spec.molecule:
                    for structure in structures:
                        if mol.isIsomorphic(structure):
                            return True, False, spec

        # At this point we can conclude that the structure does not exist
        return False, False, None

    def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True):
        """
        Formally create a new species from the specified `object`, which can be
        either a :class:`Molecule` object or an :class:`rmgpy.species.Species`
        object. It is emphasized that `reactive` relates to the :Class:`Species` attribute, while `reactive_structure`
        relates to the :Class:`Molecule` attribute.
        """

        if isinstance(object, rmgpy.species.Species):
            molecule = object.molecule[0]
            label = label if label != '' else object.label
            reactive = object.reactive
        else:
            molecule = object
            
        molecule.clearLabeledAtoms()

        # If desired, check to ensure that the species is new; return the
        # existing species if not new
        if checkForExisting:
            if isinstance(object, rmgpy.species.Species) and len(object.molecule) > 1:
                # If resonance structures were already generated (e.g., if object came from a reaction library), object
                # may contain unreactive resonance structures. Make sure a reactive structure is sent to
                # checkForExistingSpecies()
                for mol in object.molecule:
                    if mol.reactive:
                        found, reactive_structure, spec = self.checkForExistingSpecies(mol)
                        break
                else:
                    for mol in object.molecule:
                        logging.info(mol.toAdjacencyList())
                    raise AssertionError, "No reactive structures found in species {0}".format(object.molecule[0].toSMILES())
            else:
                found, reactive_structure, spec = self.checkForExistingSpecies(molecule)
            if found and reactive_structure:
                return spec, False
            if found and not reactive_structure:
                molecule.reactive=False
                spec.molecule.append(molecule)
                return spec, False

        # Check that the structure is not forbidden

        # If we're here then we're ready to make the new species
        if reactive:
            self.speciesCounter += 1   # count only reactive species
            speciesIndex = self.speciesCounter
        else:
            speciesIndex = -1
        try:
            spec = Species(index=speciesIndex, label=label, molecule=[molecule], reactive=reactive,
                 thermo=object.thermo, transportData=object.transportData)
        except AttributeError:
            spec = Species(index=speciesIndex, label=label, molecule=[molecule], reactive=reactive)
        
        spec.creationIteration = self.iterationNum
        if isinstance(object, rmgpy.species.Species) and len(object.molecule) > 1:
            # If resonance structures were already generated (e.g., if object came from a reaction library), object may
            # contain unreactive resonance structures that we'd like to keep. In this case, don't re-generate the
            # resonance structures, just keep the original ones.
            spec.molecule = object.molecule
        else:
            spec.generate_resonance_structures()
        spec.molecularWeight = Quantity(spec.molecule[0].getMolecularWeight()*1000.,"amu")
        
        if not spec.thermo:
            submit(spec,self.solventName)
        
        if spec.label == '':
            if spec.thermo and spec.thermo.label != '': #check if thermo libraries have a name for it
                logging.info('Species with SMILES of {0} named {1} based on thermo library name'.format(molecule.toSMILES().replace('/','').replace('\\',''),spec.thermo.label))
                spec.label = spec.thermo.label
                label = spec.label
            else:
                # Use SMILES as default format for label
                # However, SMILES can contain slashes (to describe the
                # stereochemistry around double bonds); since RMG doesn't 
                # distinguish cis and trans isomers, we'll just strip these out
                # so that we can use the label in file paths
                label = molecule.toSMILES().replace('/','').replace('\\','')
                
        logging.debug('Creating new species {0}'.format(label))
        
        spec.generateEnergyTransferModel()
        formula = molecule.getFormula()
        if formula in self.speciesDict:
            self.speciesDict[formula].append(spec)
        else:
            self.speciesDict[formula] = [spec]


        # Since the species is new, add it to the list of new species
        self.newSpeciesList.append(spec)

        if spec.reactive:
            self.indexSpeciesDict[spec.index] = spec

        return spec, True

    def checkForExistingReaction(self, rxn):
        """
        Check to see if an existing reaction has the same reactants, products, and
        family as `rxn`. Returns :data:`True` or :data:`False` and the matched
        reaction (if found).

        First, a shortlist of reaction is retrieved that have the same reaction keys
        as the parameter reaction.

        Next, the reaction ID containing an identifier (e.g. label) of the reactants
        and products is compared between the parameter reaction and the each of the
        reactions in the shortlist. If a match is found, the discovered reaction is 
        returned.

        If a match is not yet found, the Library (seed mechs, reaction libs)
        in the reaction database are iterated over to check if a reaction was overlooked
        (a reaction with a different "family" key as the parameter reaction).

        """

        # Make sure the reactant and product lists are sorted before performing the check
        rxn.reactants.sort()
        rxn.products.sort()

        # If reactants and products are identical, then something weird happened along
        # the way and we got a symmetrical reaction.
        if rxn.reactants == rxn.products:
            logging.debug("Symmetrical reaction found. Returning no reaction")
            return True, None
        
        familyObj = getFamilyLibraryObject(rxn.family)
        shortlist = self.searchRetrieveReactions(rxn)

        # Now use short-list to check for matches. All should be in same forward direction.

        # Make sure the reactant and product lists are sorted before performing the check
        rxn_id = generateReactionId(rxn)

        for rxn0 in shortlist:
            rxn_id0 = generateReactionId(rxn0)

            if rxn_id == rxn_id0 and areIdenticalSpeciesReferences(rxn, rxn0):
                if isinstance(familyObj, KineticsLibrary) or isinstance(familyObj, KineticsFamily):
                    if not rxn.duplicate:
                        return True, rxn0
                else:
                    return True, rxn0
            elif (isinstance(familyObj, KineticsFamily)
                  and rxn_id == rxn_id0[::-1]
                  and areIdenticalSpeciesReferences(rxn, rxn0)):
                if not rxn.duplicate:
                    return True, rxn0

        # Now check seed mechanisms
        # We want to check for duplicates in *other* seed mechanisms, but allow
        # duplicated *within* the same seed mechanism
        _, r1_fwd, r2_fwd = generateReactionKey(rxn)
        _, r1_rev, r2_rev = generateReactionKey(rxn, useProducts=True)

        for library in self.reactionDict:
            libObj = getFamilyLibraryObject(library)
            if isinstance(libObj, KineticsLibrary) and library != rxn.family:

                # First check seed short-list in forward direction                
                shortlist = self.retrieve(library, r1_fwd, r2_fwd)
                
                for rxn0 in shortlist:
                    rxn_id0 = generateReactionId(rxn0)
                    if (rxn_id == rxn_id0) or (rxn_id == rxn_id0[::-1]):
                        if areIdenticalSpeciesReferences(rxn, rxn0):
                            return True, rxn0
                
                # Now get the seed short-list of the reverse reaction

                shortlist = self.retrieve(library, r1_rev, r2_rev)
                
                for rxn0 in shortlist:
                    if areIdenticalSpeciesReferences(rxn, rxn0):
                        return True, rxn0

        return False, None

    def makeNewReaction(self, forward, checkExisting=True):
        """
        Make a new reaction given a :class:`Reaction` object `forward`. 
        The reaction is added to the global list of reactions.
        Returns the reaction in the direction that corresponds to the
        estimated kinetics, along with whether or not the reaction is new to the
        global reaction list.

        The forward direction is determined using the "is_reverse" attribute of the
        reaction's family.  If the reaction family is its own reverse, then it is
        made such that the forward reaction is exothermic at 298K.
        
        The forward reaction is appended to self.newReactionList if it is new.
        """

        # Determine the proper species objects for all reactants and products
        reactants = [self.makeNewSpecies(reactant)[0] for reactant in forward.reactants]
        products  = [self.makeNewSpecies(product)[0]  for product  in forward.products ]
        if forward.specificCollider is not None:
            forward.specificCollider = self.makeNewSpecies(forward.specificCollider)[0]

        if forward.pairs is not None:
            for pairIndex in range(len(forward.pairs)):
                reactantIndex = forward.reactants.index(forward.pairs[pairIndex][0])
                productIndex = forward.products.index(forward.pairs[pairIndex][1])
                forward.pairs[pairIndex] = (reactants[reactantIndex], products[productIndex])
                if hasattr(forward, 'reverse'):                   
                    if forward.reverse:
                        forward.reverse.pairs[pairIndex] = (products[productIndex], reactants[reactantIndex])
        forward.reactants = reactants
        forward.products  = products

        if checkExisting:
            found, rxn = self.checkForExistingReaction(forward)
            if found: return rxn, False

        # Generate the reaction pairs if not yet defined
        if forward.pairs is None:
            forward.generatePairs()
            if hasattr(forward, 'reverse'):
                if forward.reverse:
                    forward.reverse.generatePairs()
            
        # Note in the log
        if isinstance(forward, TemplateReaction):
            logging.debug('Creating new {0} template reaction {1}'.format(forward.family, forward))
        elif isinstance(forward, DepositoryReaction):
            logging.debug('Creating new {0} reaction {1}'.format(forward.getSource(), forward))
        elif isinstance(forward, LibraryReaction):
            logging.debug('Creating new library reaction {0}'.format(forward))
        else:
            raise Exception("Unrecognized reaction type {0!s}".format(forward.__class__))
        
        self.registerReaction(forward)

        forward.index = self.reactionCounter + 1
        self.reactionCounter += 1

        # Since the reaction is new, add it to the list of new reactions
        self.newReactionList.append(forward)

        # Return newly created reaction
        return forward, True

    def makeNewPDepReaction(self, forward):
        """
        Make a new pressure-dependent reaction based on a list of `reactants` and a
        list of `products`. The reaction belongs to the specified `network` and
        has pressure-dependent kinetics given by `kinetics`.

        No checking for existing reactions is made here. The returned PDepReaction
        object is not added to the global list of reactions, as that is intended
        to represent only the high-pressure-limit set. The reactionCounter is
        incremented, however, since the returned reaction can and will exist in
        the model edge and/or core.
        """

        # Don't create reverse reaction: all such reactions are treated as irreversible
        # The reverse direction will come from a different partial network
        # Note that this isn't guaranteed to satisfy thermodynamics (but will probably be close)
        forward.reverse = None
        forward.reversible = False

        # Generate the reaction pairs if not yet defined
        if forward.pairs is None:
            forward.generatePairs()
            
        # Set reaction index and increment the counter
        forward.index = self.reactionCounter + 1
        self.reactionCounter += 1

        return forward

    def enlarge(self, newObject=None, reactEdge=False,
                unimolecularReact=None, bimolecularReact=None, trimolecularReact=None):
        """
        Enlarge a reaction model by processing the objects in the list `newObject`. 
        If `newObject` is a
        :class:`rmg.species.Species` object, then the species is moved from
        the edge to the core and reactions generated for that species, reacting
        with itself and with all other species in the model core. If `newObject`
        is a :class:`rmg.unirxn.network.Network` object, then reactions are
        generated for the species in the network with the largest leak flux.

        If the `reactEdge` flag is `True`, then no newObject is needed,
        and instead the algorithm proceeds to react the core species together
        to form edge reactions.
        """
        
        numOldCoreSpecies = len(self.core.species)
        numOldCoreReactions = len(self.core.reactions)
        numOldEdgeSpecies = len(self.edge.species)
        numOldEdgeReactions = len(self.edge.reactions)
        reactionsMovedFromEdge = []
        self.newReactionList = []; self.newSpeciesList = []

        if reactEdge is False:
            # We are adding core species 
            newReactions = []
            pdepNetwork = None
            objectWasInEdge = False
        
            if isinstance(newObject, Species):
                
                newSpecies = newObject

                objectWasInEdge = newSpecies in self.edge.species
                
                if not newSpecies.reactive:
                    logging.info('NOT generating reactions for unreactive species {0}'.format(newSpecies))
                else:
                    logging.info('Adding species {0} to model core'.format(newSpecies))
                    display(newSpecies) # if running in IPython --pylab mode, draws the picture!

                # Add new species
                reactionsMovedFromEdge = self.addSpeciesToCore(newSpecies)

            elif isinstance(newObject, tuple) and isinstance(newObject[0], PDepNetwork) and self.pressureDependence:

                pdepNetwork, newSpecies = newObject
                newReactions.extend(pdepNetwork.exploreIsomer(newSpecies))

                for rxn in newReactions:
                    rxn = self.inflate(rxn)
                    try:
                        rxn.reverse = self.inflate(rxn.reverse)
                    except AttributeError:
                        pass
                    
                self.processNewReactions(newReactions, newSpecies, pdepNetwork)

            else:
                raise TypeError('Unable to use object {0} to enlarge reaction model; expecting an object of class rmg.model.Species or rmg.model.PDepNetwork, not {1}'.format(newObject, newObject.__class__))

            # If there are any core species among the unimolecular product channels
            # of any existing network, they need to be made included
            for network in self.networkList:
                network.updateConfigurations(self)
                index = 0
                isomers = [isomer.species[0] for isomer in network.isomers]
                while index < len(self.core.species):
                    species = self.core.species[index]
                    if species in isomers and species not in network.explored:
                        network.explored.append(species)
                        continue
                    for products in network.products:
                        products = products.species
                        if len(products) == 1 and products[0] == species:
                            newReactions = network.exploreIsomer(species)
                            for rxn in newReactions:
                                rxn = self.inflate(rxn)
                                try:
                                    rxn.reverse = self.inflate(rxn.reverse)
                                except AttributeError:
                                    pass

                            self.processNewReactions(newReactions, species, network)
                            network.updateConfigurations(self)
                            index = 0
                            break
                    else:
                        index += 1
            
            if isinstance(newObject, Species) and objectWasInEdge:
                # moved one species from edge to core
                numOldEdgeSpecies -= 1
                # moved these reactions from edge to core
                numOldEdgeReactions -= len(reactionsMovedFromEdge)

        else:
            # We are reacting the edge

            rxns = reactAll(self.core.species, numOldCoreSpecies,
                            unimolecularReact, bimolecularReact, trimolecularReact=trimolecularReact)
            spcs = [self.retrieveNewSpecies(rxn) for rxn in rxns]
            
            for rxn, spc in zip(rxns, spcs):
                rxn = self.inflate(rxn) 
                try:
                    rxn.reverse = self.inflate(rxn.reverse)
                except AttributeError:
                    pass
                self.processNewReactions([rxn], spc)

        ################################################################
        # Begin processing the new species and reactions
        
        # Generate kinetics of new reactions
        if self.newReactionList:
            logging.info('Generating kinetics for new reactions...')
        for reaction in self.newReactionList:
            # If the reaction already has kinetics (e.g. from a library),
            # assume the kinetics are satisfactory
            if reaction.kinetics is None:
                self.applyKineticsToReaction(reaction)
                    
        # For new reactions, convert ArrheniusEP to Arrhenius, and fix barrier heights.
        # self.newReactionList only contains *actually* new reactions, all in the forward direction.
        for reaction in self.newReactionList:
            # convert KineticsData to Arrhenius forms
            if isinstance(reaction.kinetics, KineticsData):
                reaction.kinetics = reaction.kinetics.toArrhenius()
            #  correct barrier heights of estimated kinetics
            if isinstance(reaction,TemplateReaction) or isinstance(reaction,DepositoryReaction): # i.e. not LibraryReaction
                reaction.fixBarrierHeight() # also converts ArrheniusEP to Arrhenius.
                
            if self.pressureDependence and reaction.isUnimolecular():
                # If this is going to be run through pressure dependence code,
                # we need to make sure the barrier is positive.
                reaction.fixBarrierHeight(forcePositive=True)
            
        # Update unimolecular (pressure dependent) reaction networks
        if self.pressureDependence:
            # Recalculate k(T,P) values for modified networks
            self.updateUnimolecularReactionNetworks()
            logging.info('')
            
        # Check new core and edge reactions for Chemkin duplicates
        # The same duplicate reaction gets brought into the core
        # at the same time, so there is no danger in checking all of the edge.
        newCoreReactions = self.core.reactions[numOldCoreReactions:]
        newEdgeReactions = self.edge.reactions[numOldEdgeReactions:]
        checkedReactions = self.core.reactions[:numOldCoreReactions] + self.edge.reactions[:numOldEdgeReactions]
        from rmgpy.chemkin import markDuplicateReaction
        for rxn in newCoreReactions:
            markDuplicateReaction(rxn, checkedReactions)
            checkedReactions.append(rxn)
        if self.saveEdgeSpecies:
            for rxn in newEdgeReactions:
                markDuplicateReaction(rxn, checkedReactions)
                checkedReactions.append(rxn)
        self.printEnlargeSummary(
            newCoreSpecies=self.core.species[numOldCoreSpecies:],
            newCoreReactions=self.core.reactions[numOldCoreReactions:],
            reactionsMovedFromEdge=reactionsMovedFromEdge,
            newEdgeSpecies=self.edge.species[numOldEdgeSpecies:],
            newEdgeReactions=self.edge.reactions[numOldEdgeReactions:],
            reactEdge=reactEdge,
        )

        logging.info('')
    
    def addNewSurfaceObjects(self,obj,newSurfaceSpecies,newSurfaceReactions,reactionSystem):
        """
        obj is the list of objects for enlargement coming from simulate
        newSurfaceSpecies and newSurfaceReactions are the current lists of surface species and surface reactions
        following simulation
        reactionSystem is the current reactor
        manages surface species and reactions being moved to and from the surface
        moves them to appropriate newSurfaceSpc/RxnsAdd/loss sets
        returns false if the surface has changed
        """
        surfSpcs = set(self.surface.species)
        surfRxns = set(self.surface.reactions)
        
        newSurfaceSpecies = set(newSurfaceSpecies)
        newSurfaceReactions = set(newSurfaceReactions)
                    
        addedRxns = {k for k in obj if isinstance(k,Reaction)}
        addedSurfaceRxns = newSurfaceReactions - surfRxns
                    
        addedBulkRxns = addedRxns-addedSurfaceRxns
        lostSurfaceRxns = (surfRxns - newSurfaceReactions) | addedBulkRxns
                    
        addedSpcs = {k for k in obj if isinstance(k,Species)} | {k.getMaximumLeakSpecies(reactionSystem.T.value_si, reactionSystem.P.value_si) for k in obj if isinstance(k,PDepNetwork)}
        lostSurfaceSpcs = (surfSpcs-newSurfaceSpecies) | addedSpcs
        addedSurfaceSpcs = newSurfaceSpecies - surfSpcs
                    
        self.newSurfaceSpcsAdd = self.newSurfaceSpcsAdd | addedSurfaceSpcs
        self.newSurfaceRxnsAdd = self.newSurfaceRxnsAdd | addedSurfaceRxns
        self.newSurfaceSpcsLoss = self.newSurfaceSpcsLoss | lostSurfaceSpcs
        self.newSurfaceRxnsLoss = self.newSurfaceRxnsLoss | lostSurfaceRxns
        
        return not (self.newSurfaceRxnsAdd != set() or self.newSurfaceRxnsLoss != set() or self.newSurfaceSpcsLoss != set() or self.newSurfaceSpcsAdd != set())
    
    def adjustSurface(self):
        """
        Here we add species intended to be added and remove any species that need to be moved out of the core.  
        For now we remove reactions from the surface that have become part of a PDepNetwork by 
        intersecting the set of surface reactions with the core so that all surface reactions are in the core
        thus the surface algorithm currently (June 2017) is not implemented for pdep networks
        (however it will function fine for non-pdep reactions on a pdep run)
        """
        self.surface.species = list(((set(self.surface.species) | self.newSurfaceSpcsAdd)-self.newSurfaceSpcsLoss) & set(self.core.species))
        self.surface.reactions = list(((set(self.surface.reactions) | self.newSurfaceRxnsAdd)-self.newSurfaceRxnsLoss) & set(self.core.reactions))
        self.clearSurfaceAdjustments()
        
    def clearSurfaceAdjustments(self):
        """
        empties surface tracking varaibles
        """
        self.newSurfaceSpcsAdd = set()
        self.newSurfaceRxnsAdd = set()
        self.newSurfaceSpcsLoss = set()
        self.newSurfaceRxnsLoss = set()
        
    def processNewReactions(self, newReactions, newSpecies, pdepNetwork=None):
        """
        Process a list of newly-generated reactions involving the new core
        species or explored isomer `newSpecies` in network `pdepNetwork`.
        
        Makes a reaction and decides where to put it: core, edge, or PDepNetwork.
        """
        for rxn in newReactions:
            rxn, isNew = self.makeNewReaction(rxn)
            if rxn is None:
                # Skip this reaction because there was something wrong with it
                continue
            spcs = []
            if isNew:
                # We've made a new reaction, so make sure the species involved
                # are in the core or edge
                allSpeciesInCore = True
                # Add the reactant and product species to the edge if necessary
                # At the same time, check if all reactants and products are in the core
                for spec in rxn.reactants:
                    if spec not in self.core.species:
                        allSpeciesInCore = False
                        if spec not in self.edge.species:
                            spcs.append(spec)
                            self.addSpeciesToEdge(spec)
                for spec in rxn.products:
                    if spec not in self.core.species:
                        allSpeciesInCore = False
                        if spec not in self.edge.species:
                            spcs.append(spec)
                            self.addSpeciesToEdge(spec)
                    
            isomerAtoms = sum([len(spec.molecule[0].atoms) for spec in rxn.reactants])
            
            # Decide whether or not to handle the reaction as a pressure-dependent reaction
            pdep = True
            if not self.pressureDependence:
                # The pressure dependence option is turned off entirely
                pdep = False
            elif self.pressureDependence.maximumAtoms is not None and self.pressureDependence.maximumAtoms < isomerAtoms:
                # The reaction involves so many atoms that pressure-dependent effects are assumed to be negligible
                pdep = False
            elif not (rxn.isIsomerization() or rxn.isDissociation() or rxn.isAssociation()):
                # The reaction is not unimolecular in either direction, so it cannot be pressure-dependent
                pdep = False
            elif isinstance(rxn,LibraryReaction):
                # Try generating the high pressure limit kinetics. If successful, set pdep to ``True``, and vice versa.
                pdep = rxn.generate_high_p_limit_kinetics()

            # If pressure dependence is on, we only add reactions that are not unimolecular;
            # unimolecular reactions will be added after processing the associated networks
            if not pdep:
                if not isNew: 
                    # The reaction is not new, so it should already be in the core or edge
                    continue
                if allSpeciesInCore:
                    self.addReactionToCore(rxn)
                else:
                    self.addReactionToEdge(rxn)
            else:
                # Add the reaction to the appropriate unimolecular reaction network
                # If pdepNetwork is not None then that will be the network the
                # (path) reactions are added to
                # Note that this must be done even with reactions that are not new
                # because of the way partial networks are explored
                # Since PDepReactions are created as irreversible, not doing so
                # would cause you to miss the reverse reactions!
                self.addReactionToUnimolecularNetworks(rxn, newSpecies=newSpecies, network=pdepNetwork)
                if isinstance(rxn, LibraryReaction):
                    # If reaction came from a reaction library, omit it from the core and edge so that it does 
                    # not get double-counted with the pdep network
                    if rxn in self.core.reactions:
                        self.core.reactions.remove(rxn)
                    if rxn in self.edge.reactions:
                        self.edge.reactions.remove(rxn)
            
            if not numpy.isinf(self.toleranceThermoKeepSpeciesInEdge) and spcs != []: #do thermodynamic filtering
                self.thermoFilterSpecies(spcs)
                
    def applyKineticsToReaction(self, reaction):
        """
        retrieve the best kinetics for the reaction and apply it towards the forward 
        or reverse direction (if reverse, flip the direaction).
        """
        from rmgpy.data.rmg import getDB
        # Find the reaction kinetics
        kinetics, source, entry, isForward = self.generateKinetics(reaction)
        # Flip the reaction direction if the kinetics are defined in the reverse direction
        if not isForward:
            family = getDB('kinetics').families[reaction.family]
            reaction.reactants, reaction.products = reaction.products, reaction.reactants
            reaction.pairs = [(p,r) for r,p in reaction.pairs]
            if family.ownReverse and hasattr(reaction,'reverse'):
                if reaction.reverse:
                    reaction.template = reaction.reverse.template
                    # replace degeneracy
                    reaction.degeneracy = reaction.reverse.degeneracy
                # We're done with the "reverse" attribute, so delete it to save a bit of memory
                reaction.reverse = None
        reaction.kinetics = kinetics

    def generateKinetics(self, reaction):
        """
        Generate best possible kinetics for the given `reaction` using the kinetics database.
        """
        # Only reactions from families should be missing kinetics
        assert isinstance(reaction, TemplateReaction)
        
        family = getFamilyLibraryObject(reaction.family)

        # Get the kinetics for the reaction
        kinetics, source, entry, isForward = family.getKinetics(reaction, templateLabels=reaction.template, degeneracy=reaction.degeneracy, estimator=self.kineticsEstimator, returnAllKinetics=False)
        # Get the gibbs free energy of reaction at 298 K
        G298 = reaction.getFreeEnergyOfReaction(298)
        gibbsIsPositive = G298 > -1e-8
        
        if family.ownReverse and hasattr(reaction,'reverse'):
            if reaction.reverse:
                # The kinetics family is its own reverse, so we could estimate kinetics in either direction
                
                # First get the kinetics for the other direction
                rev_kinetics, rev_source, rev_entry, rev_isForward = family.getKinetics(reaction.reverse, templateLabels=reaction.reverse.template, degeneracy=reaction.reverse.degeneracy, estimator=self.kineticsEstimator, returnAllKinetics=False)
                # Now decide which direction's kinetics to keep
                keepReverse = False
                if (entry is not None and rev_entry is None):
                    # Only the forward has an entry, meaning an exact match in a depository or template
                    # the reverse must have used an averaged estimated node - so use forward.
                    reason = "This direction matched an entry in {0}, the other was just an estimate.".format(reaction.family)
                elif (entry is None and rev_entry is not None):
                    # Only the reverse has an entry (see above) - use reverse.
                    keepReverse = True
                    reason = "This direction matched an entry in {0}, the other was just an estimate.".format(reaction.family)
                elif (entry is not None and rev_entry is not None 
                      and entry is rev_entry):
                    # Both forward and reverse have the same source and entry
                    # Use the one for which the kinetics is the forward kinetics
                    keepReverse = gibbsIsPositive and isForward and rev_isForward
                    reason = "Both directions matched the same entry in {0}, but this direction is exergonic.".format(reaction.family)
                elif self.kineticsEstimator == 'group additivity' and (kinetics.comment.find("Fitted to 1 rate")>0
                      and not rev_kinetics.comment.find("Fitted to 1 rate")>0) :
                        # forward kinetics were fitted to only 1 rate, but reverse are hopefully better
                        keepReverse = True
                        reason = "Other direction matched a group only fitted to 1 rate."
                elif self.kineticsEstimator == 'group additivity' and (not kinetics.comment.find("Fitted to 1 rate")>0
                      and rev_kinetics.comment.find("Fitted to 1 rate")>0) :
                        # reverse kinetics were fitted to only 1 rate, but forward are hopefully better
                        keepReverse = False
                        reason = "Other direction matched a group only fitted to 1 rate."
                elif entry is not None and rev_entry is not None:
                    # Both directions matched explicit rate rules
                    # Keep the direction with the lower (but nonzero) rank
                    if entry.rank < rev_entry.rank and entry.rank != 0:
                        keepReverse = False
                        reason = "Both directions matched explicit rate rules, but this direction has a rule with a lower rank ({0} vs {1}).".format(entry.rank, rev_entry.rank)
                    elif rev_entry.rank < entry.rank and rev_entry.rank != 0:
                        keepReverse = True
                        reason = "Both directions matched explicit rate rules, but this direction has a rule with a lower rank ({0} vs {1}).".format(rev_entry.rank, entry.rank)
                    # Otherwise keep the direction that is exergonic at 298 K
                    else:
                        keepReverse = gibbsIsPositive and isForward and rev_isForward
                        reason = "Both directions matched explicit rate rules, but this direction is exergonic."
                else:
                    # Keep the direction that is exergonic at 298 K
                    # This must be done after the thermo generation step
                    keepReverse = gibbsIsPositive and isForward and rev_isForward
                    reason = "Both directions are estimates, but this direction is exergonic."

                if keepReverse:
                    kinetics = rev_kinetics
                    source = rev_source
                    entry = rev_entry
                    isForward = not rev_isForward
                    G298 = -G298
                
                if self.verboseComments:
                    kinetics.comment += "\nKinetics were estimated in this direction instead of the reverse because:\n{0}".format(reason)
                    kinetics.comment += "\ndGrxn(298 K) = {0:.2f} kJ/mol".format( G298 / 1000.)
            
        # The comments generated by the database for estimated kinetics can
        # be quite long, and therefore not very useful
        # We don't want to waste lots of memory storing these long, 
        # uninformative strings, so here we replace them with much shorter ones
        if not self.verboseComments:
            # Only keep a short comment (to save memory)
            if 'Exact' in kinetics.comment:
                # Exact match of rate rule
                pass
            elif 'Matched reaction' in kinetics.comment:
                # Stems from matching a reaction from a depository
                pass
            else:
                # Estimated (averaged) rate rule
                kinetics.comment =  kinetics.comment[kinetics.comment.find('Estimated'):]
                
        return kinetics, source, entry, isForward
    
    def printEnlargeSummary(self, newCoreSpecies, newCoreReactions, newEdgeSpecies, newEdgeReactions, reactionsMovedFromEdge=None, reactEdge=False):
        """
        Output a summary of a model enlargement step to the log. The details of
        the enlargement are passed in the `newCoreSpecies`, `newCoreReactions`,
        `newEdgeSpecies`, and `newEdgeReactions` objects. 
        """

        logging.info('')
        if reactEdge:
            logging.info('Summary of Secondary Model Edge Enlargement')
        else:
            logging.info('Summary of Model Enlargement')
        logging.info('---------------------------------')

        logging.info('Added {0:d} new core species'.format(len(newCoreSpecies)))
        for spec in newCoreSpecies:
            display(spec)
            logging.info('    {0}'.format(spec))

        logging.info('Created {0:d} new edge species'.format(len(newEdgeSpecies)))
        for spec in newEdgeSpecies:
            display(spec)
            logging.info('    {0}'.format(spec))
        
        if reactionsMovedFromEdge:
            logging.info('Moved {0:d} reactions from edge to core'.format(len(reactionsMovedFromEdge)))
            for rxn in reactionsMovedFromEdge:
                for r in newCoreReactions:
                    if ((r.reactants == rxn.reactants and r.products == rxn.products) or
                        (r.products == rxn.reactants and r.reactants == rxn.products)):
                        logging.info('    {0}'.format(r))
                        newCoreReactions.remove(r)
                        break        
        logging.info('Added {0:d} new core reactions'.format(len(newCoreReactions)))
        for rxn in newCoreReactions:
            logging.info('    {0}'.format(rxn))
            
        logging.info('Created {0:d} new edge reactions'.format(len(newEdgeReactions)))
        for rxn in newEdgeReactions:
            logging.info('    {0}'.format(rxn))

        coreSpeciesCount, coreReactionCount, edgeSpeciesCount, edgeReactionCount = self.getModelSize()

        # Output current model size information after enlargement
        logging.info('')
        logging.info('After model enlargement:')
        logging.info('    The model core has {0:d} species and {1:d} reactions'.format(coreSpeciesCount, coreReactionCount))
        logging.info('    The model edge has {0:d} species and {1:d} reactions'.format(edgeSpeciesCount, edgeReactionCount))
        logging.info('')

    def addSpeciesToCore(self, spec):
        """
        Add a species `spec` to the reaction model core (and remove from edge if
        necessary). This function also moves any reactions in the edge that gain
        core status as a result of this change in status to the core.
        If this are any such reactions, they are returned in a list.
        """

        assert spec not in self.core.species, "Tried to add species {0} to core, but it's already there".format(spec.label)

        forbidden_structures = getDB('forbidden')
        
        # check RMG globally forbidden structures
        if not spec.explicitlyAllowed and forbidden_structures.isMoleculeForbidden(spec.molecule[0]):
            
            rxnList = []
            if spec in self.edge.species:

                #remove forbidden species from edge
                logging.info("Species {0} was Forbidden and not added to Core...Removing from Edge.".format(spec))
                self.removeSpeciesFromEdge(self.reactionSystems,spec)

                return []
        
        # Add the species to the core
        self.core.species.append(spec)
        
        rxnList = []
        if spec in self.edge.species:

            # If species was in edge, remove it
            logging.debug("Removing species {0} from edge.".format(spec))
            self.edge.species.remove(spec)

            # Search edge for reactions that now contain only core species;
            # these belong in the model core and will be moved there
            for rxn in self.edge.reactions:
                allCore = True
                for reactant in rxn.reactants:
                    if reactant not in self.core.species: allCore = False
                for product in rxn.products:
                    if product not in self.core.species: allCore = False
                if allCore: rxnList.append(rxn)

            # Move any identified reactions to the core
            for rxn in rxnList:
                self.addReactionToCore(rxn)
                logging.debug("Moving reaction from edge to core: {0}".format(rxn))
        return rxnList

    def addSpeciesToEdge(self, spec):
        """
        Add a species `spec` to the reaction model edge.
        """
        self.edge.species.append(spec)
    
    def setThermodynamicFilteringParameters(self,Tmax, toleranceThermoKeepSpeciesInEdge,minCoreSizeForPrune,maximumEdgeSpecies,reactionSystems):
        """
        sets parameters for thermodynamic filtering based on the current core
        Tmax is the maximum reactor temperature in K
        toleranceThermoKeepSpeciesInEdge is the Gibbs number above which species will be filtered
        minCoreSizeForPrune is the core size at which thermodynamic filtering will start
        maximumEdgeSpecies is the maximum allowed number of edge species
        reactionSystems is a list of reactionSystem objects
        """
        self.Tmax = Tmax
        Gs = [spc.thermo.getFreeEnergy(Tmax) for spc in self.core.species]
        self.Gmax = max(Gs)
        self.Gmin = min(Gs)
        
        self.Gfmax = toleranceThermoKeepSpeciesInEdge*(self.Gmax-self.Gmin)+self.Gmax
        self.toleranceThermoKeepSpeciesInEdge = toleranceThermoKeepSpeciesInEdge
        self.minCoreSizeForPrune = minCoreSizeForPrune
        self.reactionSystems = reactionSystems
        self.maximumEdgeSpecies = maximumEdgeSpecies
    
    def thermoFilterSpecies(self, spcs):
        """
        checks Gibbs energy of the species in species against the
        maximum allowed Gibbs energy
        """
        Tmax = self.Tmax
        for spc in spcs:
            G = spc.thermo.getFreeEnergy(Tmax)
            if G > self.Gfmax:
                Gn = (G-self.Gmax)/(self.Gmax-self.Gmin)
                logging.info('Removing species {0} with Gibbs energy {1} from edge because it\'s Gibbs number {2} is greater than the toleranceThermoKeepSpeciesInEdge of {3} '.format(spc,G,Gn,self.toleranceThermoKeepSpeciesInEdge))
                self.removeSpeciesFromEdge(self.reactionSystems,spc)
                
        # Delete any networks that became empty as a result of pruning
        if self.pressureDependence:
            self.removeEmptyPdepNetworks()
                        
    def thermoFilterDown(self,maximumEdgeSpecies,minSpeciesExistIterationsForPrune=0):
        """
        removes species from the edge based on their Gibbs energy until maximumEdgeSpecies
        is reached under the constraint that all removed species are older than
        minSpeciesExistIterationsForPrune iterations
        maximumEdgeSpecies is the maximum allowed number of edge species
        minSpeciesExistIterationsForPrune is the number of iterations a species must be in the edge
        before it is eligible for thermo filtering
        """
        Tmax = self.Tmax
        numToRemove = len(self.edge.species) - maximumEdgeSpecies
        logging.debug('Planning to remove {0} species'.format(numToRemove))
        iteration = self.iterationNum
            
        if numToRemove > 0: #implies flux pruning is off or did not trigger
            logging.info('Reached maximum number of edge species')
            logging.info('Attempting to remove excess edge species with Thermodynamic filtering')
            spcs = self.edge.species
            Gfs = numpy.array([spc.thermo.getFreeEnergy(Tmax) for spc in spcs])
            Gns = (Gfs-self.Gmax)/(self.Gmax-self.Gmin) 
            inds = numpy.argsort(Gns) #could actually do this with the Gfs, but want to print the Gn value later
            inds = inds[::-1] #get in order of increasing Gf

            ind = 0
            removeSpcs = []
            
            rInds = []
            while ind < len(inds) and numToRemove > 0: #find the species we can remove and collect indices for removal     
                i = inds[ind]
                spc = spcs[i]
                if iteration - spc.creationIteration >= minSpeciesExistIterationsForPrune:
                    removeSpcs.append(spc)
                    rInds.append(i)
                    numToRemove -= 1
                ind += 1
            
            logging.debug('found {0} eligible species for filtering'.format(len(removeSpcs)))
            
            for i,spc in enumerate(removeSpcs):
                logging.info('Removing species {0} from edge to meet maximum number of edge species, Gibbs number is {1}'.format(spc,Gns[rInds[i]]))
                self.removeSpeciesFromEdge(self.reactionSystems,spc)
            
            # Delete any networks that became empty as a result of pruning
            if self.pressureDependence:
                self.removeEmptyPdepNetworks()
            
            #call garbage collection
            collected = gc.collect()
            logging.info('Garbage collector: collected %d objects.' % (collected))
            
    def removeEmptyPdepNetworks(self):
        """
        searches for and deletes any empty pdep networks
        """
        networksToDelete = []
        for network in self.networkList:
            if len(network.pathReactions) == 0 and len(network.netReactions) == 0:
                networksToDelete.append(network)
                    
        if len(networksToDelete) > 0:
            logging.info('Deleting {0:d} empty pressure-dependent reaction networks'.format(len(networksToDelete)))
            for network in networksToDelete:
                logging.debug('    Deleting empty pressure dependent reaction network #{0:d}'.format(network.index))
                source = tuple(network.source)
                nets_with_this_source = self.networkDict[source]
                nets_with_this_source.remove(network)
                if not nets_with_this_source:
                    del(self.networkDict[source])
                self.networkList.remove(network)
                    
    def prune(self, reactionSystems, toleranceKeepInEdge, toleranceMoveToCore, maximumEdgeSpecies, minSpeciesExistIterationsForPrune):
        """
        Remove species from the model edge based on the simulation results from
        the list of `reactionSystems`.
        """

        ineligibleSpecies = []     # A list of the species which are not eligible for pruning, for any reason
        prunableSpecies = reactionSystems[0].prunableSpecies
        prunableNetworks = reactionSystems[0].prunableNetworks
        
        numPrunableSpecies = len(prunableSpecies)
        iteration = self.iterationNum
        # All edge species that have not existed for more than two enlarge
        # iterations are ineligible for pruning
        for spec in prunableSpecies:
            if iteration - spec.creationIteration <= minSpeciesExistIterationsForPrune:
                ineligibleSpecies.append(spec)

        # Get the maximum species rates (and network leak rates)
        # across all reaction systems
        maxEdgeSpeciesRateRatios = numpy.zeros((numPrunableSpecies), numpy.float64)
        for reactionSystem in reactionSystems:
            for i in range(numPrunableSpecies):
                rateRatio = reactionSystem.maxEdgeSpeciesRateRatios[i]
                if maxEdgeSpeciesRateRatios[i] < rateRatio:
                    maxEdgeSpeciesRateRatios[i] = rateRatio

            for i,network in enumerate(prunableNetworks):
                rateRatio = reactionSystem.maxNetworkLeakRateRatios[i]
                # Add the fraction of the network leak rate contributed by
                # each unexplored species to that species' rate
                # This is to ensure we have an overestimate of that species flux
                ratios = network.getLeakBranchingRatios(reactionSystem.T.value_si,reactionSystem.P.value_si)
                for spec, frac in ratios.iteritems():
                    if spec in prunableSpecies:
                        index = prunableSpecies.index(spec)
                        maxEdgeSpeciesRateRatios[index] += frac * rateRatio
                # Mark any species that is explored in any partial network as ineligible for pruning
                for spec in network.explored:
                    if spec not in ineligibleSpecies:
                        ineligibleSpecies.append(spec)

        # Sort the edge species rates by index
        indices = numpy.argsort(maxEdgeSpeciesRateRatios)
        # Determine which species to prune
        speciesToPrune = []
        pruneDueToRateCounter = 0
        for index in indices:
            spec = prunableSpecies[index]
            if spec in ineligibleSpecies or not spec in self.edge.species:
                continue
            # Remove the species with rates below the pruning tolerance from the model edge
            if maxEdgeSpeciesRateRatios[index] < toleranceKeepInEdge:
                speciesToPrune.append((index, spec))
                pruneDueToRateCounter += 1
            # Keep removing species with the lowest rates until we are below the maximum edge species size
            elif numPrunableSpecies - len(speciesToPrune) > maximumEdgeSpecies:
                if maxEdgeSpeciesRateRatios[index] < toleranceMoveToCore:
                    logging.info('Pruning species {0} to make numEdgeSpecies smaller than maximumEdgeSpecies'.format(spec))
                    speciesToPrune.append((index, spec))
                else:
                    logging.warning('Attempted to prune a species that exceeded toleranceMoveToCore, pruning settings for this run are likely bad, either maximumEdgeSpecies needs to be set higher (~100000) or minSpeciesExistIterationsForPrune should be reduced (~2)')
                    break
            else:
                break

        # Actually do the pruning
        if pruneDueToRateCounter > 0:
            logging.info('Pruning {0:d} species whose rate ratios against characteristic rate did not exceed the minimum threshold of {1:g}'.format(pruneDueToRateCounter, toleranceKeepInEdge))
            for index, spec in speciesToPrune[0:pruneDueToRateCounter]:
                logging.info('Pruning species {0:<56}'.format(spec))
                logging.debug('    {0:<56}    {1:10.4e}'.format(spec, maxEdgeSpeciesRateRatios[index]))
                self.removeSpeciesFromEdge(reactionSystems, spec)
        if len(speciesToPrune) - pruneDueToRateCounter > 0:
            logging.info('Pruning {0:d} species to obtain an edge size of {1:d} species'.format(len(speciesToPrune) - pruneDueToRateCounter, maximumEdgeSpecies))
            for index, spec in speciesToPrune[pruneDueToRateCounter:]:
                logging.info('Pruning species {0:<56}'.format(spec))
                logging.debug('    {0:<56}    {1:10.4e}'.format(spec, maxEdgeSpeciesRateRatios[index]))
                self.removeSpeciesFromEdge(reactionSystems, spec)

        # Delete any networks that became empty as a result of pruning
        if self.pressureDependence:
            self.removeEmptyPdepNetworks()

        logging.info('')


    def removeSpeciesFromEdge(self, reactionSystems, spec):
        """
        Remove species `spec` from the reaction model edge.
        """

        # remove the species
        self.edge.species.remove(spec)
        self.indexSpeciesDict.pop(spec.index)

        # clean up species references in reactionSystems
        for reactionSystem in reactionSystems:
            try:
                reactionSystem.speciesIndex.pop(spec)
            except KeyError:
                pass

            # identify any reactions it's involved in
            rxnList = []
            for rxn in reactionSystem.reactionIndex:
                if spec in rxn.reactants or spec in rxn.products:
                    rxnList.append(rxn)

            for rxn in rxnList:
                reactionSystem.reactionIndex.pop(rxn)

        # identify any reactions it's involved in
        rxnList = []
        for rxn in self.edge.reactions:
            if spec in rxn.reactants or spec in rxn.products:
                rxnList.append(rxn)
        # remove those reactions
        for rxn in rxnList:
            self.edge.reactions.remove(rxn)
        
        # Remove the species from any unirxn networks it is in
        if self.pressureDependence:
            for network in self.networkList:
                # Delete all path reactions involving the species
                rxnList = []
                for rxn in network.pathReactions:
                    if spec in rxn.reactants or spec in rxn.products:
                        rxnList.append(rxn)
                if len(rxnList) > 0:
                    for rxn in rxnList:
                        network.pathReactions.remove(rxn)
                    # Delete all net reactions involving the species
                    rxnList = []
                    for rxn in network.netReactions:
                        if spec in rxn.reactants or spec in rxn.products:
                            rxnList.append(rxn)
                    for rxn in rxnList:
                        network.netReactions.remove(rxn)
                        
                    # Recompute the isomers, reactants, and products for this network
                    network.updateConfigurations(self)

        # Remove from the global list of reactions
        # also remove it from the global list of reactions
        for family in self.reactionDict:
            if spec in self.reactionDict[family]:
                del self.reactionDict[family][spec]
            for reactant1 in self.reactionDict[family]:
                if spec in self.reactionDict[family][reactant1]:
                    del self.reactionDict[family][reactant1][spec]
            for reactant1 in self.reactionDict[family]:
                for reactant2 in self.reactionDict[family][reactant1]:
                    tempRxnDeleteList = []
                    for templateReaction in self.reactionDict[family][reactant1][reactant2]:
                        if spec in templateReaction.reactants or spec in templateReaction.products:
                            tempRxnDeleteList.append(templateReaction)
                    for tempRxnToBeDeleted in tempRxnDeleteList:
                        self.reactionDict[family][reactant1][reactant2].remove(tempRxnToBeDeleted)

        # remove from the global list of species, to free memory
        formula = spec.molecule[0].getFormula()
        self.speciesDict[formula].remove(spec)
        if spec in self.speciesCache:
            self.speciesCache.remove(spec)
            self.speciesCache.append(None)

    def addReactionToCore(self, rxn):
        """
        Add a reaction `rxn` to the reaction model core (and remove from edge if
        necessary). This function assumes `rxn` has already been checked to
        ensure it is supposed to be a core reaction (i.e. all of its reactants
        AND all of its products are in the list of core species).
        """
        if rxn not in self.core.reactions:
            self.core.reactions.append(rxn)
        if rxn in self.edge.reactions:
            self.edge.reactions.remove(rxn)
        
    def addReactionToEdge(self, rxn):
        """
        Add a reaction `rxn` to the reaction model edge. This function assumes
        `rxn` has already been checked to ensure it is supposed to be an edge
        reaction (i.e. all of its reactants OR all of its products are in the
        list of core species, and the others are in either the core or the
        edge).
        """
        self.edge.reactions.append(rxn)

    def getModelSize(self):
        """
        Return the numbers of species and reactions in the model core and edge.
        Note that this is not necessarily equal to the lengths of the
        corresponding species and reaction lists.
        """
        coreSpeciesCount = len(self.core.species)
        coreReactionsCount = len(self.core.reactions)
        edgeSpeciesCount = len(self.edge.species)
        edgeReactionsCount = len(self.edge.reactions)
        return (coreSpeciesCount, coreReactionsCount, edgeSpeciesCount, edgeReactionsCount)

    def getLists(self):
        """
        Return lists of all of the species and reactions in the core and the
        edge.
        """
        speciesList = []
        speciesList.extend(self.core.species)
        speciesList.extend(self.edge.species)
        reactionList = []
        reactionList.extend(self.core.reactions)
        reactionList.extend(self.edge.reactions)
        return speciesList, reactionList

    def getStoichiometryMatrix(self):
        """
        Return the stoichiometry matrix for all generated species and reactions.
        The id of each species and reaction is the corresponding row and column,
        respectively, in the matrix.
        """
        speciesList, reactionList = self.getLists()
        from scipy import sparse
        stoichiometry = sparse.dok_matrix((self.speciesCounter, self.reactionCounter), float)
        for rxn in reactionList:
            j = rxn.index - 1
            specList = rxn.reactants[:]; specList.extend(rxn.products)
            for spec in specList:
                i = spec.index - 1
                nu = rxn.getStoichiometricCoefficient(spec)
                if nu != 0: stoichiometry[i,j] = nu
        return stoichiometry.tocsr()

    def addSeedMechanismToCore(self, seedMechanism, react=False):
        """
        Add all species and reactions from `seedMechanism`, a 
        :class:`KineticsPrimaryDatabase` object, to the model core. If `react`
        is ``True``, then reactions will also be generated between the seed
        species. For large seed mechanisms this can be prohibitively expensive,
        so it is not done by default.
        """

        if react: raise NotImplementedError("react=True doesn't work yet")
        database = rmgpy.data.rmg.database

        libraryNames = database.kinetics.libraries.keys()
        familyNames = database.kinetics.families.keys()
        
        path = os.path.join(settings['database.directory'],'kinetics','libraries')
        from rmgpy.rmg.input import rmg
        
        self.newReactionList = []; self.newSpeciesList = []

        numOldCoreSpecies = len(self.core.species)
        numOldCoreReactions = len(self.core.reactions)

        logging.info('Adding seed mechanism {0} to model core...'.format(seedMechanism))

        seedMechanism = database.kinetics.libraries[seedMechanism]
        
        rxns = seedMechanism.getLibraryReactions()
        
        for rxn in rxns:
            if isinstance(rxn,LibraryReaction) and not (rxn.library in libraryNames) and not (rxn.library == 'kineticsjobs'): #if one of the reactions in the library is from another library load that library
                database.kinetics.libraryOrder.append((rxn.library,'Internal'))
                database.kinetics.loadLibraries(path=path,libraries=[rxn.library])
                libraryNames = database.kinetics.libraries.keys()
            if isinstance(rxn,TemplateReaction) and not (rxn.family in familyNames):
                logging.warning('loading reaction {0} originally from family {1} as a library reaction'.format(str(rxn),rxn.family))
                rxn = LibraryReaction(reactants=rxn.reactants[:], products=rxn.products[:],
                 library=seedMechanism.name, specificCollider=rxn.specificCollider, kinetics=rxn.kinetics, duplicate=rxn.duplicate,
                 reversible=rxn.reversible
                 )
            r, isNew = self.makeNewReaction(rxn) # updates self.newSpeciesList and self.newReactionlist
            if not isNew:
                logging.info("This library reaction was not new: {0}".format(rxn))
            elif self.pressureDependence and rxn.elementary_high_p and rxn.isUnimolecular()\
                    and isinstance(rxn, LibraryReaction) and isinstance(rxn.kinetics, Arrhenius):
                # This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics.
                # We should calculate a pressure-dependent rate for it
                if len(rxn.reactants) == 1:
                    self.processNewReactions(newReactions=[rxn],newSpecies=rxn.reactants[0])
                else:
                    self.processNewReactions(newReactions=[rxn],newSpecies=rxn.products[0])
            
        # Perform species constraints and forbidden species checks
        
        for spec in self.newSpeciesList:
            if database.forbiddenStructures.isMoleculeForbidden(spec.molecule[0]):
                if 'allowed' in rmg.speciesConstraints and 'seed mechanisms' in rmg.speciesConstraints['allowed']:
                    spec.explicitlyAllowed = True
                    logging.warning("Species {0} from seed mechanism {1} is globally forbidden.  It will behave as an inert unless found in a seed mechanism or reaction library.".format(spec.label, seedMechanism.label))
                else:
                    raise ForbiddenStructureException("Species {0} from seed mechanism {1} is globally forbidden. You may explicitly allow it, but it will remain inert unless found in a seed mechanism or reaction library.".format(spec.label, seedMechanism.label))
            if failsSpeciesConstraints(spec):
                if 'allowed' in rmg.speciesConstraints and 'seed mechanisms' in rmg.speciesConstraints['allowed']:
                    rmg.speciesConstraints['explicitlyAllowedMolecules'].extend(spec.molecule)
                else:
                    raise ForbiddenStructureException("Species constraints forbids species {0} from seed mechanism {1}. Please reformulate constraints, remove the species, or explicitly allow it.".format(spec.label, seedMechanism.label))

        for spec in self.newSpeciesList:            
            if spec.reactive:
                submit(spec,self.solventName)

            self.addSpeciesToCore(spec)

        for rxn in self.newReactionList:
            if self.pressureDependence and rxn.isUnimolecular():
                # If this is going to be run through pressure dependence code,
                # we need to make sure the barrier is positive.
                # ...but are Seed Mechanisms run through PDep? Perhaps not.
                for spec in itertools.chain(rxn.reactants, rxn.products):
                    submit(spec,self.solventName)

                rxn.fixBarrierHeight(forcePositive=True)
            self.addReactionToCore(rxn)
        
        # Check we didn't introduce unmarked duplicates
        self.markChemkinDuplicates()
        
        self.printEnlargeSummary(
            newCoreSpecies=self.core.species[numOldCoreSpecies:],
            newCoreReactions=self.core.reactions[numOldCoreReactions:],
            newEdgeSpecies=[],
            newEdgeReactions=[],
        )



    def addReactionLibraryToEdge(self, reactionLibrary):
        """
        Add all species and reactions from `reactionLibrary`, a
        :class:`KineticsPrimaryDatabase` object, to the model edge.
        """

        database = rmgpy.data.rmg.database
        libraryNames = database.kinetics.libraries.keys()
        familyNames = database.kinetics.families.keys()
        path = os.path.join(settings['database.directory'],'kinetics','libraries')
        
        from rmgpy.rmg.input import rmg

        self.newReactionList = []
        self.newSpeciesList = []

        numOldEdgeSpecies = len(self.edge.species)
        numOldEdgeReactions = len(self.edge.reactions)

        logging.info('Adding reaction library {0} to model edge...'.format(reactionLibrary))
        reactionLibrary = database.kinetics.libraries[reactionLibrary]

        rxns = reactionLibrary.getLibraryReactions()
        for rxn in rxns:
            if isinstance(rxn,LibraryReaction) and not (rxn.library in libraryNames): #if one of the reactions in the library is from another library load that library
                database.kinetics.libraryOrder.append((rxn.library,'Internal'))
                database.kinetics.loadLibraries(path=path,libraries=[rxn.library])
                libraryNames = database.kinetics.libraries.keys()
            if isinstance(rxn,TemplateReaction) and not (rxn.family in familyNames):
                logging.warning('loading reaction {0} originally from family {1} as a library reaction'.format(str(rxn),rxn.family))
                rxn = LibraryReaction(reactants=rxn.reactants[:], products=rxn.products[:],
                 library=reactionLibrary.name, specificCollider=rxn.specificCollider, kinetics=rxn.kinetics, duplicate=rxn.duplicate,
                 reversible=rxn.reversible
                 )
            r, isNew = self.makeNewReaction(rxn) # updates self.newSpeciesList and self.newReactionlist
            if not isNew:
                logging.info("This library reaction was not new: {0}".format(rxn))
            elif self.pressureDependence and rxn.elementary_high_p and rxn.isUnimolecular()\
                    and isinstance(rxn, LibraryReaction) and isinstance(rxn.kinetics, Arrhenius):
                # This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics.
                # We should calculate a pressure-dependent rate for it
                if len(rxn.reactants) == 1:
                    self.processNewReactions(newReactions=[rxn],newSpecies=rxn.reactants[0])
                else:
                    self.processNewReactions(newReactions=[rxn],newSpecies=rxn.products[0])

        # Perform species constraints and forbidden species checks
        for spec in self.newSpeciesList:
            if database.forbiddenStructures.isMoleculeForbidden(spec.molecule[0]):
                if 'allowed' in rmg.speciesConstraints and 'reaction libraries' in rmg.speciesConstraints['allowed']:
                    spec.explicitlyAllowed = True
                    logging.warning("Species {0} from reaction library {1} is globally forbidden.  It will behave as an inert unless found in a seed mechanism or reaction library.".format(spec.label, reactionLibrary.label))
                else:
                    raise ForbiddenStructureException("Species {0} from reaction library {1} is globally forbidden. You may explicitly allow it, but it will remain inert unless found in a seed mechanism or reaction library.".format(spec.label, reactionLibrary.label))
            if failsSpeciesConstraints(spec):
                if 'allowed' in rmg.speciesConstraints and 'reaction libraries' in rmg.speciesConstraints['allowed']:
                    rmg.speciesConstraints['explicitlyAllowedMolecules'].extend(spec.molecule)
                else:
                    raise ForbiddenStructureException("Species constraints forbids species {0} from reaction library {1}. Please reformulate constraints, remove the species, or explicitly allow it.".format(spec.label, reactionLibrary.label))

        for spec in self.newSpeciesList:
            if spec.reactive: 
                submit(spec,self.solventName)

            self.addSpeciesToEdge(spec)

        for rxn in self.newReactionList:
            # Note that we haven't actually evaluated any fluxes at this point
            # Instead, we remove the comment below if the reaction is moved to
            # the core later in the mechanism generation
            if not (self.pressureDependence and rxn.elementary_high_p and rxn.isUnimolecular()
                    and isinstance(rxn, LibraryReaction) and isinstance(rxn.kinetics, Arrhenius)):
                # Don't add to the edge library reactions that were already processed
                self.addReactionToEdge(rxn)

        if self.saveEdgeSpecies:
            from rmgpy.chemkin import markDuplicateReaction
            newEdgeReactions = self.edge.reactions[numOldEdgeReactions:]
            checkedReactions = self.core.reactions + self.edge.reactions[:numOldEdgeReactions]
            for rxn in newEdgeReactions:
                markDuplicateReaction(rxn, checkedReactions)
                checkedReactions.append(rxn)

        self.printEnlargeSummary(
            newCoreSpecies=[],
            newCoreReactions=[],
            newEdgeSpecies=self.edge.species[numOldEdgeSpecies:],
            newEdgeReactions=self.edge.reactions[numOldEdgeReactions:],
        )

    def addReactionLibraryToOutput(self, reactionLib):
        """
        Add all species and reactions from `reactionLibrary`, a
        :class:`KineticsPrimaryDatabase` object, to the output.
        This does not bring any of the reactions or species into the core itself.
        """

        logging.info('Adding reaction library {0} to output file...'.format(reactionLib))
        
        # Append the edge reactions that are from the selected reaction library to an output species and output reactions list
        for rxn in self.edge.reactions:
            if isinstance(rxn, LibraryReaction):
                if rxn.library == reactionLib:
                    self.outputReactionList.append(rxn)
                    
                    for species in rxn.reactants + rxn.products:
                        if species not in self.core.species and species not in self.outputSpeciesList:
                            self.outputSpeciesList.append(species)
                            


    def addReactionToUnimolecularNetworks(self, newReaction, newSpecies, network=None):
        """
        Given a newly-created :class:`Reaction` object `newReaction`, update the
        corresponding unimolecular reaction network. If no network exists, a new
        one is created. If the new reaction is an isomerization that connects two
        existing networks, the two networks are merged. This function is called
        whenever a new high-pressure limit edge reaction is created. Returns the
        network containing the new reaction.
        """

        assert isinstance(newSpecies, Species)

        # Put the reaction in the direction in which the new species is in the reactants
        if newSpecies in newReaction.reactants:
            reactants = newReaction.reactants[:]
            products = newReaction.products[:]
        else:
            reactants = newReaction.products[:]
            products = newReaction.reactants[:] 
        reactants.sort()
        products.sort()
        
        source = tuple(reactants)

        # Only search for a network if we don't specify it as a parameter
        if network is None:
            if len(reactants) == 1:
                # Find the network containing the reactant as the source
                try:
                    networks = self.networkDict[source]
                    assert len(networks) == 1
                    network = networks[0]
                except KeyError:
                    pass
            elif len(reactants) > 1:
                # Find the network containing the reactants as the source AND the
                # product as an explored isomer
                try:
                    networks = self.networkDict[source]
                    for n in networks:
                        if products[0] in n.explored:
                            assert network is None
                            network = n
                except KeyError:
                    pass
            else:
                return

            # If no suitable network exists, create a new one
            if network is None:
                self.networkCount += 1
                network = PDepNetwork(index=self.networkCount, source=reactants[:])
                        # should the source passed to PDepNetwork constuctor be a tuple not a list? that's what is used in networkDict
                try:
                    self.networkDict[source].append(network)
                except KeyError:
                    self.networkDict[source] = [network]
                self.networkList.append(network)

        # Add the path reaction to that network
        network.addPathReaction(newReaction)

    def updateUnimolecularReactionNetworks(self):
        """
        Iterate through all of the currently-existing unimolecular reaction
        networks, updating those that have been marked as invalid. In each update,
        the phenomonological rate coefficients :math:`k(T,P)` are computed for
        each net reaction in the network, and the resulting reactions added or
        updated.
        """

        # Merge networks if necessary
        # Two partial networks having the same source and containing one or
        # more explored isomers in common must be merged together to avoid
        # double-counting of rates
        for networks in self.networkDict.itervalues():
            networkCount = len(networks)
            for index0, network0 in enumerate(networks):
                index = index0 + 1
                while index < networkCount:
                    found = False
                    network = networks[index]
                    if network0.source == network.source:
                        # The networks contain the same source, but do they contain any common included isomers (other than the source)?
                        for isomer in network0.explored:
                            if isomer != network.source and isomer in network.explored:
                                # The networks contain an included isomer in common, so we need to merge them
                                found = True
                                break
                    if found:
                        # The networks contain the same source and one or more common included isomers
                        # Therefore they need to be merged together
                        logging.info('Merging PDepNetwork #{0:d} and PDepNetwork #{1:d}'.format(network0.index, network.index))
                        network0.merge(network)
                        networks.remove(network)
                        self.networkList.remove(network)
                        networkCount -= 1
                    else:
                        index += 1

        count = sum([1 for network in self.networkList if not network.valid and not (len(network.explored) == 0 and len(network.source) > 1)])
        logging.info('Updating {0:d} modified unimolecular reaction networks (out of {1:d})...'.format(count, len(self.networkList)))
        
        # Iterate over all the networks, updating the invalid ones as necessary
        # self = reactionModel object
        updatedNetworks = []
        for network in self.networkList:
            if not network.valid:
                network.update(self, self.pressureDependence)
                updatedNetworks.append(network)
            
        # PDepReaction objects generated from partial networks are irreversible
        # However, it makes more sense to have reversible reactions in the core
        # Thus we mark PDepReaction objects as reversible and remove the reverse
        # direction from the list of core reactions
        # Note that well-skipping reactions may not have a reverse if the well
        # that they skip over is not itself in the core
        index = 0
        coreReactionCount = len(self.core.reactions)
        while index < coreReactionCount:
            reaction = self.core.reactions[index]
            if isinstance(reaction, PDepReaction):
                for reaction2 in self.core.reactions[index+1:]:
                    if isinstance(reaction2, PDepReaction) and reaction.reactants == reaction2.products and reaction.products == reaction2.reactants:
                        # We've found the PDepReaction for the reverse direction
                        dGrxn = reaction.getFreeEnergyOfReaction(300.)
                        kf = reaction.getRateCoefficient(1000,1e5)
                        kr = reaction.getRateCoefficient(1000,1e5) / reaction.getEquilibriumConstant(1000)
                        kf2 = reaction2.getRateCoefficient(1000,1e5) / reaction2.getEquilibriumConstant(1000)
                        kr2 = reaction2.getRateCoefficient(1000,1e5)
                        if kf / kf2 < 0.5 or kf / kf2 > 2.0:
                            # Most pairs of reactions should satisfy thermodynamic consistency (or at least be "close")
                            # Warn about the ones that aren't close (but don't abort)
                            logging.warning('Forward and reverse PDepReactions for reaction {0!s} generated from networks {1:d} and {2:d} do not satisfy thermodynamic consistency.'.format(reaction, reaction.network.index, reaction2.network.index))
                            logging.warning('{0!s}:'.format(reaction))
                            logging.warning('{0:.2e} {1:.2e}:'.format(kf, kf2))
                            logging.warning('{0!s}:'.format(reaction2))
                            logging.warning('{0:.2e} {1:.2e}:'.format(kr, kr2))
                        # Keep the exergonic direction
                        keepFirst = dGrxn < 0
                        # Delete the PDepReaction that we aren't keeping
                        if keepFirst:
                            self.core.reactions.remove(reaction2)
                            reaction.reversible = True
                        else:
                            self.core.reactions.remove(reaction)
                            self.core.reactions.remove(reaction2)
                            self.core.reactions.insert(index, reaction2)
                            reaction2.reversible = True
                        coreReactionCount -= 1
                        # There should be only one reverse, so we can stop searching once we've found it
                        break
                else:
                    reaction.reversible = True
            # Move to the next core reaction
            index += 1


    def markChemkinDuplicates(self):
        """
        Check that all reactions that will appear the chemkin output have been checked as duplicates.
        
        Call this if you've done something that may have introduced undetected duplicate reactions,
        like add a reaction library or seed mechanism.
        Anything added via the :meth:`expand` method should already be detected.
        """
        from rmgpy.chemkin import markDuplicateReactions
        
        rxnList = self.core.reactions + self.outputReactionList
        markDuplicateReactions(rxnList)
        
    
    def registerReaction(self, rxn):
        """
        Adds the reaction to the reaction database.

        The reaction database is structured as a multi-level
        dictionary, for efficient search and retrieval of
        existing reactions.

        The database has two types of dictionary keys:
        - reaction family
        - reactant(s) keys

        First, the keys are generated for the parameter reaction.
        
        Next, it is checked whether the reaction database already 
        contains similar keys. If not, a new container is created,
        either a dictionary for the family key and first reactant key,
        or a list for the second reactant key.

        Finally, the reaction is inserted as the first element in the 
        list.
        """

        key_family, key1, key2 = generateReactionKey(rxn)

        # make dictionary entries if necessary
        if key_family not in self.reactionDict:
            self.reactionDict[key_family] = {}

        if not self.reactionDict[key_family].has_key(key1):
            self.reactionDict[key_family][key1] = {}

        if not self.reactionDict[key_family][key1].has_key(key2):
            self.reactionDict[key_family][key1][key2] = []

        # store this reaction at the top of the relevant short-list
        self.reactionDict[key_family][key1][key2].insert(0, rxn)


    def searchRetrieveReactions(self, rxn):
        """
        Searches through the reaction database for 
        reactions with an identical reaction key as the key of the 
        parameter reaction.

        Both the reaction key based on the reactants as well as on the products
        is used to search for possible candidate reactions.
        """

        # Get the short-list of reactions with the same family, reactant1 and reactant2
        family_label, r1_fwd, r2_fwd = generateReactionKey(rxn)
        
        my_reactionList = []

        rxns = self.retrieve(family_label, r1_fwd, r2_fwd)
        my_reactionList.extend(rxns)
            
            
        family = getFamilyLibraryObject(family_label)       
        # if the family is its own reverse (H-Abstraction) then check the other direction
        if isinstance(family,KineticsFamily): 

            # Get the short-list of reactions with the same family, product1 and product2
            family_label, r1_rev, r2_rev = generateReactionKey(rxn, useProducts=True)

            rxns = self.retrieve(family_label, r1_rev, r2_rev)
            my_reactionList.extend(rxns)

        return my_reactionList

    def initializeIndexSpeciesDict(self):
        """
        Populates the core species dictionary
        
        integer -> core Species

        with the species that are currently in the core.
        """

        for spc in itertools.chain(self.core.species, self.edge.species):
            if spc.reactive:
                self.indexSpeciesDict[spc.index] = spc

    def retrieve(self, family_label, key1, key2):
        """
        Returns a list of reactions from the reaction database with the 
        same keys as the parameters.

        Returns an empty list when one of the keys could not be found.
        """
        try:
            return self.reactionDict[family_label][key1][key2][:]
        except KeyError: # no such short-list: must be new, unless in seed.
            return []

    def inflate(self, rxn):
        """
        Convert reactions from
        reactants/products that are referring
        to the core species index, to the respective Species objects.
        """
        reactants, products, pairs = [], [], []

        for reactant in rxn.reactants:
            reactant = self.getSpecies(reactant)  
            reactants.append(reactant)

        for product in rxn.products:
            product = self.getSpecies(product)
            products.append(product)

        for reactant, product in rxn.pairs:
            reactant = self.getSpecies(reactant)
            product = self.getSpecies(product)
            pairs.append((reactant, product))

        rxn.reactants = reactants  
        rxn.products = products
        rxn.pairs = pairs

        return rxn

    def getSpecies(self, obj):
        """
        Retrieve species object, by
        polling the index species dictionary.
        """
        if isinstance(obj, int):
            spc = self.indexSpeciesDict[obj]
            return spc
        return obj

    def retrieveNewSpecies(self, deflatedRxn):
        """
        Searches for the first reactant or product in the deflated reaction
        that is represented by an integer.

        Such an object refers to a core species that was used to generate the
        reaction in the first place. Reactants or products represented by an
        object that is not an integer will be a newly-generated structure.
        """
        for obj in itertools.chain(deflatedRxn.reactants, deflatedRxn.products):
            if isinstance(obj, int):
                return self.getSpecies(obj)
        raise Exception("No core species were found in either reactants or products of {0}!".format(deflatedRxn))


def generateReactionKey(rxn, useProducts=False):
    """
    Returns a tuple with 3 keys:
    - the reaction family (or library) the reaction belongs to
    - the keys of the reactants.

    None for the third element in the tuple if there is 
    only 1 reactant.

    The keys are sorted alphabetically.
    """

    key_family = rxn.family

    spc_list = rxn.products if useProducts else rxn.reactants
    key1 = getKey(spc_list[0])
    key2 = None if len(spc_list) == 1 else getKey(spc_list[1])
    key1, key2 = sorted([key1, key2], reverse=True)# ensure None is always at end

    return (key_family, key1, key2)

def generateReactionId(rxn):
    """
    Returns a tuple of the reactions reactant and product
    keys.

    Both lists are sorted.

    The first element in the tuple is the reactants list.
    """


    reactants = sorted([getKey(reactant) for reactant in rxn.reactants])
    products = sorted([getKey(product) for product in rxn.products])

    return (reactants, products)

def getFamilyLibraryObject(label):
    """
    Returns the KineticsFamily or KineticsLibrary object associated with the
    parameter string.

    First search through the reaction families, then 
    through the libraries.
    """

    kinetics = rmgpy.data.rmg.database.kinetics

    try:
        fam = kinetics.families[label]
        return fam
    except KeyError:
        pass

    try:
        lib = kinetics.libraries[label]
        return lib
    except KeyError:
        pass

    raise Exception('Could not retrieve the family/library: {}'.format(label))

def getKey(spc):
    """
    Returns a string of the species that can serve as a key in a dictionary.
    """

    return spc.label

def areIdenticalSpeciesReferences(rxn1, rxn2):
    """
    Checks if the references of the reactants and products of the two reactions
    are identical, in either direction.
    """
    identical_same_direction = rxn1.reactants == rxn2.reactants and rxn1.products == rxn2.products
    identical_opposite_directions = rxn1.reactants == rxn2.products and rxn1.products == rxn2.reactants
    identical_collider = rxn1.specificCollider == rxn2.specificCollider
    
    return (identical_same_direction or identical_opposite_directions) and identical_collider
